{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PyQuante in Julia\n",
    "Experimenting with writing quantum chemistry in Julia"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "trace2 (generic function with 1 method)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "factorial2(n::Int64) = prod(n:-2:1) # double factorial !!\n",
    "dist2(dx,dy,dz) = dx*dx+dy*dy+dz*dz # Is there something in the standard library that does this?\n",
    "\n",
    "pairs(n::Int64) = ((i, j) for i = 1:n for j = 1:i)\n",
    "rpairs(n::Int64) = ((i,j) for i in 1:n for j in 1:n) # rectangular option to old pairs\n",
    "spairs(n::Int64) = ((i, j) for i = 1:n for j = 1:(i-1)) # subdiagonal option to old pairs\n",
    " \n",
    "triangle(i::Int64) = div(i*(i+1),2)\n",
    "triangle(i::Int64,j::Int64) = i<j ? triangle(j-1)+i : triangle(i-1)+j\n",
    "                        \n",
    "iiterator(n::Int64) = ((i,j,k,l) for (i,j) in pairs(n) for (k,l) in pairs(n) if triangle(i,j) <= triangle(k,l))\n",
    "\n",
    "iindex(i::Int64,j::Int64,k::Int64,l::Int64) = triangle(triangle(i,j),triangle(k,l))\n",
    "trace2(A,B) = sum(A.*B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "@assert factorial2(6)==48\n",
    "@assert collect(pairs(3)) == [(1,1),(2,1),(2,2),(3,1),(3,2),(3,3)]\n",
    "#@assert collect(pairs(3,\"subdiag\")) == [(2,1),(3,1),(3,2)]\n",
    "@assert collect(spairs(3))== [(2,1),(3,1),(3,2)]\n",
    "#@assert collect(pairs(2,\"rect\")) == [(1,1),(1,2),(2,1),(2,2)]\n",
    "@assert collect(rpairs(2)) == [(1,1),(1,2),(2,1),(2,2)]\n",
    "@assert iindex(1,1,1,1) == 1\n",
    "@assert iindex(1,1,1,2) == iindex(1,1,2,1) == iindex(1,2,1,1) == iindex(2,1,1,1) == 2\n",
    "@assert iindex(1,1,2,2) == iindex(2,2,1,1) == 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basis function definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "normalize! (generic function with 1 method)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type PGBF\n",
    "    expn::Float64\n",
    "    x::Float64\n",
    "    y::Float64\n",
    "    z::Float64\n",
    "    I::Int64\n",
    "    J::Int64\n",
    "    K::Int64\n",
    "    norm::Float64\n",
    "end\n",
    "\n",
    "function pgbf(expn,x=0,y=0,z=0,I=0,J=0,K=0,norm=1)\n",
    "    p = PGBF(expn,x,y,z,I,J,K,norm)\n",
    "    normalize!(p)\n",
    "    return p\n",
    "end\n",
    "\n",
    "function amplitude(bf::PGBF,x,y,z)\n",
    "    dx,dy,dz = x-bf.x,y-bf.y,z-bf.z\n",
    "    r2 = dist2(dx,dy,dz)\n",
    "    return bf.norm*(dx^bf.I)*(dy^bf.J)*(dz^bf.K)*exp(-bf.expn*r2)\n",
    "end\n",
    "\n",
    "function normalize!(pbf::PGBF)\n",
    "    pbf.norm /= sqrt(overlap(pbf,pbf))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "contract (generic function with 2 methods)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type CGBF\n",
    "    x::Float64\n",
    "    y::Float64\n",
    "    z::Float64\n",
    "    I::Int64\n",
    "    J::Int64\n",
    "    K::Int64\n",
    "    norm::Float64\n",
    "    pgbfs::Array{PGBF,1}\n",
    "    coefs::Array{Float64,1}\n",
    "end\n",
    "\n",
    "cgbf(x=0,y=0,z=0,I=0,J=0,K=0) = CGBF(x,y,z,I,J,K,1.0,PGBF[],Float64[])\n",
    "\n",
    "function amplitude(bf::CGBF,x,y,z)\n",
    "    s = 0\n",
    "    for (c,pbf) in primitives(bf)\n",
    "        s += c*amplitude(pbf,x,y,z)\n",
    "    end\n",
    "    return bf.norm*s\n",
    "end\n",
    "\n",
    "function normalize!(bf::CGBF)\n",
    "    bf.norm /= sqrt(overlap(bf,bf))\n",
    "end\n",
    "\n",
    "primitives(a::CGBF) = zip(a.coefs,a.pgbfs)\n",
    "\n",
    "function push!(cbf::CGBF,expn,coef)\n",
    "    Base.push!(cbf.pgbfs,pgbf(expn,cbf.x,cbf.y,cbf.z,cbf.I,cbf.J,cbf.K))\n",
    "    Base.push!(cbf.coefs,coef)\n",
    "    normalize!(cbf)\n",
    "end\n",
    "\n",
    "function contract(f,a::CGBF,b::CGBF)\n",
    "    s = 0\n",
    "    for (ca,abf) in primitives(a)\n",
    "        for (cb,bbf) in primitives(b)\n",
    "            s += ca*cb*f(abf,bbf)\n",
    "        end\n",
    "    end\n",
    "    return a.norm*b.norm*s\n",
    "end\n",
    "\n",
    "function contract(f,a::CGBF,b::CGBF,c::CGBF,d::CGBF)\n",
    "    s = 0\n",
    "    for (ca,abf) in primitives(a)\n",
    "        for (cb,bbf) in primitives(b)\n",
    "            for (cc,cbf) in primitives(c)\n",
    "                for (cd,dbf) in primitives(d)\n",
    "                    s += ca*cb*cc*cd*f(abf,bbf,cbf,dbf)\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    return a.norm*b.norm*c.norm*d.norm*s\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = pgbf(1.0)\n",
    "px = pgbf(1.0,0,0,0,1,0,0)\n",
    "@assert isapprox(amplitude(s,0,0,0),0.71270547)\n",
    "@assert isapprox(amplitude(px,0,0,0),0)\n",
    "c = cgbf(0.0,0.0,0.0)\n",
    "push!(c,1,1)\n",
    "@assert isapprox(amplitude(c,0,0,0),0.71270547)\n",
    "c2 = cgbf(0,0,0)\n",
    "push!(c2,1,0.2)\n",
    "push!(c2,0.5,0.2)\n",
    "@assert isapprox(overlap(c2,c2),1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## One-electron integrals\n",
    "### Overlap matrix elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "binomial_prefactor (generic function with 2 methods)"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function overlap(a::PGBF,b::PGBF)\n",
    "    return a.norm*b.norm*overlap(a.expn,a.x,a.y,a.z,a.I,a.J,a.K,\n",
    "    b.expn,b.x,b.y,b.z,b.I,b.J,b.K)\n",
    "end\n",
    "\n",
    "overlap(a::CGBF,b::CGBF) = contract(overlap,a,b)\n",
    "\n",
    "function overlap(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK)\n",
    "    gamma = aexpn+bexpn\n",
    "    px,py,pz = gaussian_product_center(aexpn,ax,ay,az,bexpn,bx,by,bz)\n",
    "    rab2 = dist2(ax-bx,ay-by,az-bz) \n",
    "    pre = (pi/gamma)^1.5*exp(-aexpn*bexpn*rab2/gamma)\n",
    "    wx = overlap1d(aI,bI,px-ax,px-bx,gamma)\n",
    "    wy = overlap1d(aJ,bJ,py-ay,py-by,gamma)\n",
    "    wz = overlap1d(aK,bK,pz-az,pz-bz,gamma)\n",
    "    return pre*wx*wy*wz\n",
    "end\n",
    "\n",
    "function gaussian_product_center(a::PGBF,b::PGBF)\n",
    "    return (a.expn*[a.x,a.y,a.z]+b.expn*[b.x,b.y,b.z])/(a.expn+b.expn)\n",
    "end\n",
    "\n",
    "function gaussian_product_center(aexpn::Float64,ax::Float64,ay::Float64,az::Float64,\n",
    "        bexpn::Float64,bx::Float64,by::Float64,bz::Float64)\n",
    "    return (aexpn*[ax,ay,az]+bexpn*[bx,by,bz])/(aexpn+bexpn)    \n",
    "end\n",
    "\n",
    "function overlap1d(la::Int64,lb::Int64,ax::Float64,bx::Float64,gamma::Float64)\n",
    "    total = 0\n",
    "    for i in 0:div(la+lb,2)\n",
    "        total += binomial_prefactor(2i,la,lb,ax,bx)*factorial2(2i-1)/(2gamma)^i\n",
    "    end\n",
    "    return total\n",
    "end\n",
    "\n",
    "function binomial_prefactor(s::Int64,ia::Int64,ib::Int64,xpa::Float64,xpb::Float64)\n",
    "    #println(\"binomial_prefactor($s,$ia,$ib,$xpa,$xpb)\")\n",
    "    total = 0\n",
    "    for t in 0:s\n",
    "        if (s-ia) <= t <= ib\n",
    "            total += binomial(ia,s-t)*binomial(ib,t)*xpa^(ia-s+t)*xpb^(ib-t)\n",
    "        end\n",
    "    end\n",
    "    return total\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "@assert overlap1d(0,0,0,0,1) == 1\n",
    "@assert gaussian_product_center(s,s) == [0,0,0]\n",
    "@assert isapprox(overlap(s,s),1)\n",
    "@assert isapprox(overlap(px,px),1)\n",
    "@assert isapprox(overlap(s,px),0)\n",
    "@assert binomial_prefactor(0,0,0,0,0) == 1\n",
    "@assert isapprox(overlap(c,c),1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kinetic matrix elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "kinetic (generic function with 3 methods)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function kinetic(a::PGBF,b::PGBF)\n",
    "    return a.norm*b.norm*kinetic(a.expn,a.x,a.y,a.z,a.I,a.J,a.K,\n",
    "                                b.expn,b.x,b.y,b.z,b.I,b.J,b.K)\n",
    "end\n",
    "\n",
    "function kinetic(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK)\n",
    "    overlap0 = overlap(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK)\n",
    "    overlapx1 = overlap(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI+2,bJ,bK)\n",
    "    overlapy1 = overlap(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ+2,bK)\n",
    "    overlapz1 = overlap(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK+2)\n",
    "    overlapx2 = overlap(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI-2,bJ,bK)\n",
    "    overlapy2 = overlap(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ-2,bK)\n",
    "    overlapz2 = overlap(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK-2)\n",
    "    term0 = bexpn*(2*(bI+bJ+bK)+3)*overlap0\n",
    "    term1 = -2*(bexpn^2)*(overlapx1+overlapy1+overlapz1)\n",
    "    term2 = -0.5*(bI*(bI-1)*overlapx2+bJ*(bJ-1)*overlapy2+bK*(bK-1)*overlapz2)\n",
    "    return term0+term1+term2\n",
    "end\n",
    "\n",
    "kinetic(a::CGBF,b::CGBF) = contract(kinetic,a,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "@assert isapprox(kinetic(1,0,0,0,0,0,0,1,0,0,0,0,0,0),2.9530518648229536)\n",
    "@assert isapprox(kinetic(s,s),1.5)\n",
    "@assert isapprox(kinetic(c,c),1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nuclear attraction term"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "nuclear_attraction (generic function with 9 methods)"
      ]
     },
     "execution_count": 167,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function Aterm(i,r,u,l1,l2,ax,bx,cx,gamma)\n",
    "    term1 = (-1)^i*binomial_prefactor(i,l1,l2,ax,bx)\n",
    "    term2 = (-1)^u*factorial(i)*cx^(i-2r-2u)\n",
    "    term3 = (1/4/gamma)^(r+u)/factorial(r)/factorial(u)/factorial(i-2r-2u)\n",
    "    return term1*term2*term3\n",
    "end\n",
    "\n",
    "function Aarray(l1,l2,a,b,c,g)\n",
    "    Imax = l1+l2+1\n",
    "    A = zeros(Float64,Imax)\n",
    "    for i in 0:(Imax-1)\n",
    "        for r in 0:div(i,2)\n",
    "            for u in 0:div(i-2r,2)\n",
    "                I = i-2r-u+1\n",
    "                A[I] += Aterm(i,r,u,l1,l2,a,b,c,g)\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    return A\n",
    "end\n",
    "\n",
    "function nuclear_attraction(aexpn::Float64,ax::Float64,ay::Float64,az::Float64,\n",
    "                            aI::Int64,aJ::Int64,aK::Int64,\n",
    "                            bexpn::Float64,bx::Float64,by::Float64,bz::Float64,\n",
    "                            bI::Int64,bJ::Int64,bK::Int64,\n",
    "                            cx::Float64,cy::Float64,cz::Float64)\n",
    "    px,py,pz = gaussian_product_center(aexpn,ax,ay,az,bexpn,bx,by,bz)\n",
    "    gamma = aexpn+bexpn\n",
    "    rab2 = dist2(ax-bx,ay-by,az-bz)\n",
    "    rcp2 = dist2(cx-px,cy-py,cz-pz)\n",
    "    Ax = Aarray(aI,bI,px-ax,px-bx,px-cx,gamma)\n",
    "    Ay = Aarray(aJ,bJ,py-ay,py-by,py-cy,gamma)\n",
    "    Az = Aarray(aK,bK,pz-az,pz-bz,pz-cz,gamma)\n",
    "    total = 0\n",
    "    for I in 0:(aI+bI)\n",
    "        for J in 0:(aJ+bJ)\n",
    "            for K in 0:(aK+bK)\n",
    "                total += Ax[I+1]*Ay[J+1]*Az[K+1]*Fgamma(I+J+K,rcp2*gamma)\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    val=-2pi*exp(-aexpn*bexpn*rab2/gamma)*total/gamma\n",
    "    #println(val)\n",
    "    #println((Ax,Ay,Az,rcp2*gamma,Fgamma(0,rcp2*gamma)))\n",
    "    return val\n",
    "end\n",
    "\n",
    "function nuclear_attraction(a::PGBF,b::PGBF,cx::Float64,cy::Float64,cz::Float64)\n",
    "    return a.norm*b.norm*nuclear_attraction(a.expn,a.x,a.y,a.z,a.I,a.J,a.K,\n",
    "                                            b.expn,b.x,b.y,b.z,b.I,b.J,b.K,cx,cy,cz)\n",
    "end\n",
    "nuclear_attraction(a::PGBF,b::PGBF,c::Atom) = c.atno*nuclear_attraction(a,b,c.x,c.y,c.z)\n",
    "nuclear_attraction(a::PGBF,b::PGBF,m::Molecule) = sum([nuclear_attraction(a,b,c) for c in m.atomlist])\n",
    "\n",
    "function Fgamma(m,x,SMALL=1e-12)\n",
    "    #println(\"Fgamma($m,$x)\")\n",
    "    x = max(x,SMALL) # Evidently needs underflow protection\n",
    "    return 0.5*x^(-m-0.5)*gammainc(m+0.5,x)\n",
    "end\n",
    "\n",
    "function gammainc(a::Float64,x::Float64)\n",
    "    # This is the series version of gamma from pyquante. For reasons I don't get, it \n",
    "    # doesn't work around a=1. This works alright, but is only a stopgap solution\n",
    "    # until Julia gets an incomplete gamma function programmed\n",
    "    if abs(a-1) < 1e-3\n",
    "        println(\"Warning: gammainc_series is known to have problems for a ~ 1\")\n",
    "    end\n",
    "    if x < (a+1.0)\n",
    "        #Use the series representation\n",
    "        gam,gln = gser(a,x)\n",
    "    else \n",
    "        #Use continued fractions\n",
    "        gamc,gln = gcf(a,x)\n",
    "        gam = 1-gamc\n",
    "    end\n",
    "    return exp(gln)*gam\n",
    "end\n",
    "\n",
    "function gser(a,x,ITMAX=100,EPS=3e-9)\n",
    "    # Series representation of Gamma. NumRec sect 6.1.\n",
    "    gln=lgamma(a)\n",
    "    if x == 0\n",
    "        return 0,gln\n",
    "    end\n",
    "    ap = a\n",
    "    delt = s = 1/a\n",
    "    for i in 1:ITMAX\n",
    "        ap += 1\n",
    "        delt *= (x/ap)\n",
    "        s += delt\n",
    "        if abs(delt) < abs(s)*EPS\n",
    "            break\n",
    "        end\n",
    "    end\n",
    "    return s*exp(-x+a*log(x)-gln),gln\n",
    "end\n",
    "\n",
    "function gcf(a::Float64,x::Float64,ITMAX::Int64=200,EPS::Float64=3e-9,FPMIN::Float64=1e-30)\n",
    "    #Continued fraction representation of Gamma. NumRec sect 6.1\"\n",
    "    gln=lgamma(a)\n",
    "    b=x+1.-a\n",
    "    c=1./FPMIN\n",
    "    d=1./b\n",
    "    h=d\n",
    "    for i in 1:ITMAX\n",
    "        an=-i*(i-a)\n",
    "        b=b+2.\n",
    "        d=an*d+b\n",
    "        if abs(d) < FPMIN\n",
    "            d=FPMIN\n",
    "        end\n",
    "        c=b+an/c\n",
    "        if abs(c) < FPMIN\n",
    "            c=FPMIN\n",
    "        end\n",
    "        d=1./d\n",
    "        delt=d*c\n",
    "        h=h*delt\n",
    "        if abs(delt-1.) < EPS\n",
    "            break\n",
    "        end\n",
    "    end\n",
    "    gammcf = exp(-x+a*log(x)-gln)*h\n",
    "    return gammcf,gln\n",
    "end\n",
    "\n",
    "# Need a nested scope to squeeze this into the contract function\n",
    "function nuclear_attraction(a::CGBF,b::CGBF,cx,cy,cz)\n",
    "    na(a,b) = nuclear_attraction(a,b,cx,cy,cz)\n",
    "    contract(na,a,b)\n",
    "end\n",
    "function nuclear_attraction(a::CGBF,b::CGBF,c::Atom)\n",
    "    na(a,b) = nuclear_attraction(a,b,c)\n",
    "    contract(na,a,b)\n",
    "end\n",
    "function nuclear_attraction(a::CGBF,b::CGBF,m::Molecule)\n",
    "    na(a,b) = nuclear_attraction(a,b,m)\n",
    "    contract(na,a,b)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [],
   "source": [
    "@assert isapprox(nuclear_attraction(s,s,0,0,0),-1.59576912)\n",
    "@assert isapprox(nuclear_attraction(c,c,0,0,0),-1.59576912)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [],
   "source": [
    "li,h = lih.atomlist\n",
    "bfs = build_basis(lih)\n",
    "s1,s2,x,y,z,h1s = bfs.bfs\n",
    "@assert isapprox(nuclear_attraction(s1,s1,lih),-8.307532656)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {},
   "outputs": [],
   "source": [
    "@assert Aterm(0,0,0,0,0,0,0,0,0) == 1.0\n",
    "@assert Aarray(0,0,0,0,0,1) == [1.0]\n",
    "@assert Aarray(0,1,1,1,1,1) == [1.0, -1.0]\n",
    "@assert Aarray(1,1,1,1,1,1) == [1.5, -2.5, 1.0]\n",
    "@assert Aterm(0,0,0,0,0,0,0,0,1) == 1.0\n",
    "@assert Aterm(0,0,0,0,1,1,1,1,1) == 1.0\n",
    "@assert Aterm(1,0,0,0,1,1,1,1,1) == -1.0\n",
    "@assert Aterm(0,0,0,1,1,1,1,1,1) == 1.0\n",
    "@assert Aterm(1,0,0,1,1,1,1,1,1) == -2.0\n",
    "@assert Aterm(2,0,0,1,1,1,1,1,1) == 1.0\n",
    "@assert Aterm(2,0,1,1,1,1,1,1,1) == -0.5\n",
    "@assert Aterm(2,1,0,1,1,1,1,1,1) == 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "# gammainc test functions. Test values taken from Mathematica\n",
    "# println(\"a=0.5 test\")\n",
    "@assert maximum([gammainc(0.5,x) for x in 0:10]\n",
    "        -[0, 1.49365, 1.69181, 1.7471, 1.76416, 1.76968, \n",
    "            1.77151, 1.77213, 1.77234, 1.77241, 1.77244]) < 1e-5\n",
    "    \n",
    "# println(\"a=1.5 test\")\n",
    "@assert maximum([gammainc(1.5,x) for x in 0:10]\n",
    "        -[0, 1.49365, 1.69181, 1.7471, 1.76416, 1.76968, \n",
    "            1.77151, 1.77213, 1.77234, 1.77241, 1.77244]) < 1e-5\n",
    "# println(\"a=2.5 test\")\n",
    "@assert maximum([gammainc(2.5,x) for x in 0:10]\n",
    "        -[0, 0.200538, 0.59898, 0.922271, 1.12165, 1.22933, \n",
    "            1.2831, 1.30859, 1.32024, 1.32542, 1.32768]) < 1e-5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two electron integrals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "coulomb (generic function with 3 methods)"
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function coulomb(aexpn,ax,ay,az,aI,aJ,aK,\n",
    "    bexpn,bx,by,bz,bI,bJ,bK,\n",
    "    cexpn,cx,cy,cz,cI,cJ,cK,\n",
    "    dexpn,dx,dy,dz,dI,dJ,dK)\n",
    "    # This is the slow method of computing integrals from Huzinaga et al.\n",
    "    # Use the HRR/VRR scheme from Head-Gordon & Pople instead\n",
    "\n",
    "    rab2 = dist2(ax-bx,ay-by,az-bz)\n",
    "    rcd2 = dist2(cx-dx,cy-dy,cz-dz)\n",
    "    \n",
    "    px,py,pz = gaussian_product_center(aexpn,ax,ay,az,bexpn,bx,by,bz)\n",
    "    qx,qy,qz = gaussian_product_center(cexpn,cx,cy,cz,dexpn,dx,dy,dz)\n",
    "    rpq2 = dist2(px-qx,py-qy,pz-qz)\n",
    "\n",
    "    g1 = aexpn+bexpn\n",
    "    g2 = cexpn+dexpn\n",
    "    delta = 0.25*(1/g1+1/g2)\n",
    "    \n",
    "    Bx = Barray(aI,bI,cI,dI,px,ax,bx,qx,cx,dx,g1,g2,delta)\n",
    "    By = Barray(aJ,bJ,cJ,dJ,py,ay,by,qy,cy,dy,g1,g2,delta)\n",
    "    Bz = Barray(aK,bK,cK,dK,pz,az,bz,qz,cz,dz,g1,g2,delta)\n",
    "    \n",
    "    s = 0\n",
    "    for I in 0:(aI+bI+cI+dI)\n",
    "        for J in 0:(aJ+bJ+cJ+dJ)\n",
    "            for K in 0:(aK+bK+cK+dK)\n",
    "                s += Bx[I+1]*By[J+1]*Bz[K+1]*Fgamma(I+J+K,0.25*rpq2/delta)\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    return 2*pi^(2.5)/(g1*g2*sqrt(g1+g2))*exp(-aexpn*bexpn*rab2/g1)*exp(-cexpn*dexpn*rcd2/g2)*s\n",
    "end\n",
    "\n",
    "function coulomb(a::PGBF,b::PGBF,c::PGBF,d::PGBF)\n",
    "    return a.norm*b.norm*c.norm*d.norm*coulomb(a.expn,a.x,a.y,a.z,a.I,a.J,a.K,\n",
    "        b.expn,b.x,b.y,b.z,b.I,b.J,b.K,\n",
    "        c.expn,c.x,c.y,c.z,c.I,c.J,c.K,\n",
    "        d.expn,d.x,d.y,d.z,d.I,d.J,d.K)\n",
    "end\n",
    "\n",
    "fB(i::Int64,l1::Int64,l2::Int64,p::Float64,a::Float64,b::Float64,r::Float64,g::Float64) = binomial_prefactor(i,l1,l2,p-a,p-b)*B0(i,r,g)\n",
    "B0(i::Int64,r::Int64,g::Float64) = fact_ratio2(i,r)*(4g)^(r-i)\n",
    "fact_ratio2(a::Int64,b::Int64) = factorial(a)/factorial(b)/factorial(a-2b)\n",
    "\n",
    "function Bterm(i1::Int64,i2::Int64,r1::Int64,r2::Int64,u::Int64,l1::Int64,l2::Int64,l3::Int64,l4::Int64,\n",
    "        Px::Float64,Ax::Float64,Bx::Float64,Qx::Float64,Cx::Float64,Dx::Float64,\n",
    "        gamma1::Float64,gamma2::Float64,delta::Float64)\n",
    "    # THO eq. 2.22\n",
    "    return fB(i1,l1,l2,Px,Ax,Bx,r1,gamma1)*(-1)^i2*fB(i2,l3,l4,Qx,Cx,Dx,r2,gamma2)\n",
    "           *(-1)^u*fact_ratio2(i1+i2-2*(r1+r2),u)\n",
    "           *(Qx-Px)^(i1+i2-2*(r1+r2)-2*u)/delta^(i1+i2-2*(r1+r2)-u)\n",
    "end\n",
    "\n",
    "function Barray(l1::Int64,l2::Int64,l3::Int64,l4::Int64,p::Float64,a::Float64,b::Float64,q::Float64,\n",
    "        c::Float64,d::Float64,g1::Float64,g2::Float64,delta::Float64)\n",
    "    Imax = l1+l2+l3+l4+1\n",
    "    B = zeros(Float64,Imax)\n",
    "    for i1 in 0:(l1+l2)\n",
    "        for i2 in 0:(l3+l4)\n",
    "            for r1 in 0:div(i1,2)\n",
    "                for r2 in 0:div(i2,2)\n",
    "                    for u in 0:(div(i1+i2,2)-r1-r2)\n",
    "                        I = i1+i2-2*(r1+r2)-u\n",
    "                        B[I+1] += Bterm(i1,i2,r1,r2,u,l1,l2,l3,l4,p,a,b,q,c,d,g1,g2,delta)\n",
    "                    end\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    return B\n",
    "end\n",
    "\n",
    "coulomb(a::CGBF,b::CGBF,c::CGBF,d::CGBF) = contract(coulomb,a,b,c,d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [],
   "source": [
    "@assert fB(0,0,0,0.0,0.0,0.0,0,2.0) == 1\n",
    "@assert fB(0,0,0,1.0,1.0,1.0,0,2.0) == 1\n",
    "@assert B0(0,0,2.0) == 1\n",
    "@assert fact_ratio2(0,0) == 1\n",
    "@assert Bterm(0,0,0,0,0,0,0,0,0,0.0,0.0,0.0,0.0,0.0,0.0,2.0,2.0,0.25)==1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [],
   "source": [
    "@assert isapprox(coulomb(1, 0,0,0, 0,0,0, 1, 0,0,0, 0,0,0, 1, 0,0,0, 0,0,0, 1, 0,0,0, 0,0,0),4.37335458)\n",
    "@assert isapprox(coulomb(s,s,s,s),1.128379167)\n",
    "@assert isapprox(coulomb(c,c,c,c),1.128379167)\n",
    "@assert isapprox(coulomb(1, 0,0,0, 0,0,0, 1, 0,0,1, 0,0,0, 1, 0,0,0, 0,0,0, 1, 0,0,1, 0,0,0),1.6088672396)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basis Set Data\n",
    "Note use of curly braces here. Julia assumes that if you have square braces, you want\n",
    "things flattened as much as possible (to be as fast as possible, I guess). Curlys \n",
    "preserve the list structure the way I would expect from Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dict{String,Array{Array{Tuple{Char,Array{Tuple{Float64,Float64},1}},1},1}} with 1 entry:\n",
       "  \"sto3g\" => Array{Tuple{Char,Array{Tuple{Float64,Float64},1}},1}[Tuple{Char,Ar…"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sto3g = [\n",
    "    # H\n",
    "    [('S',\n",
    "      [(3.4252509099999999, 0.15432897000000001),\n",
    "       (0.62391373000000006, 0.53532813999999995),\n",
    "       (0.16885539999999999, 0.44463454000000002)])],\n",
    "    # He\n",
    "    [('S',\n",
    "      [(6.3624213899999997, 0.15432897000000001),\n",
    "       (1.1589229999999999, 0.53532813999999995),\n",
    "       (0.31364978999999998, 0.44463454000000002)])],\n",
    "    # Li\n",
    "    [('S',\n",
    "      [(16.119575000000001, 0.15432897000000001),\n",
    "       (2.9362007000000001, 0.53532813999999995),\n",
    "       (0.79465050000000004, 0.44463454000000002)]),\n",
    "     ('S',\n",
    "      [(0.63628969999999996, -0.099967230000000004),\n",
    "       (0.14786009999999999, 0.39951282999999999),\n",
    "       (0.048088699999999998, 0.70011546999999996)]),\n",
    "     ('P',\n",
    "      [(0.63628969999999996, 0.15591627),\n",
    "       (0.14786009999999999, 0.60768372000000004),\n",
    "       (0.048088699999999998, 0.39195739000000002)])],\n",
    "    # Be\n",
    "    [('S',\n",
    "      [(30.167871000000002, 0.15432897000000001),\n",
    "       (5.4951153000000001, 0.53532813999999995),\n",
    "       (1.4871927, 0.44463454000000002)]),\n",
    "     ('S',\n",
    "      [(1.3148331, -0.099967230000000004),\n",
    "       (0.3055389, 0.39951282999999999),\n",
    "       (0.099370700000000006, 0.70011546999999996)]),\n",
    "     ('P',\n",
    "      [(1.3148331, 0.15591627),\n",
    "       (0.3055389, 0.60768372000000004),\n",
    "       (0.099370700000000006, 0.39195739000000002)])],\n",
    "    # B\n",
    "    [('S',\n",
    "      [(48.791113000000003, 0.15432897000000001),\n",
    "       (8.8873622000000001, 0.53532813999999995),\n",
    "       (2.4052669999999998, 0.44463454000000002)]),\n",
    "     ('S',\n",
    "      [(2.2369561, -0.099967230000000004),\n",
    "       (0.51982050000000002, 0.39951282999999999),\n",
    "       (0.16906180000000001, 0.70011546999999996)]),\n",
    "     ('P',\n",
    "      [(2.2369561, 0.15591627),\n",
    "       (0.51982050000000002, 0.60768372000000004),\n",
    "       (0.16906180000000001, 0.39195739000000002)])],\n",
    "    # C\n",
    "    [('S',\n",
    "      [(71.616837000000004, 0.15432897000000001),\n",
    "       (13.045095999999999, 0.53532813999999995),\n",
    "       (3.5305122, 0.44463454000000002)]),\n",
    "     ('S',\n",
    "      [(2.9412493999999998, -0.099967230000000004),\n",
    "       (0.68348310000000001, 0.39951282999999999),\n",
    "       (0.22228990000000001, 0.70011546999999996)]),\n",
    "     ('P',\n",
    "      [(2.9412493999999998, 0.15591627),\n",
    "       (0.68348310000000001, 0.60768372000000004),\n",
    "       (0.22228990000000001, 0.39195739000000002)])],\n",
    "    # N\n",
    "    [('S',\n",
    "      [(99.106168999999994, 0.15432897000000001),\n",
    "       (18.052312000000001, 0.53532813999999995),\n",
    "       (4.8856602000000002, 0.44463454000000002)]),\n",
    "     ('S',\n",
    "      [(3.7804559000000002, -0.099967230000000004),\n",
    "       (0.87849659999999996, 0.39951282999999999),\n",
    "       (0.28571439999999998, 0.70011546999999996)]),\n",
    "     ('P',\n",
    "      [(3.7804559000000002, 0.15591627),\n",
    "       (0.87849659999999996, 0.60768372000000004),\n",
    "       (0.28571439999999998, 0.39195739000000002)])],\n",
    "    # O\n",
    "    [('S',\n",
    "      [(130.70931999999999, 0.15432897000000001),\n",
    "       (23.808861, 0.53532813999999995),\n",
    "       (6.4436083000000002, 0.44463454000000002)]),\n",
    "     ('S',\n",
    "      [(5.0331513000000001, -0.099967230000000004),\n",
    "       (1.1695960999999999, 0.39951282999999999),\n",
    "       (0.38038899999999998, 0.70011546999999996)]),\n",
    "     ('P',\n",
    "      [(5.0331513000000001, 0.15591627),\n",
    "       (1.1695960999999999, 0.60768372000000004),\n",
    "       (0.38038899999999998, 0.39195739000000002)])],\n",
    "    # F\n",
    "    [('S',\n",
    "      [(166.67912999999999, 0.15432897000000001),\n",
    "       (30.360811999999999, 0.53532813999999995),\n",
    "       (8.2168206999999995, 0.44463454000000002)]),\n",
    "     ('S',\n",
    "      [(6.4648032000000004, -0.099967230000000004),\n",
    "       (1.5022812000000001, 0.39951282999999999),\n",
    "       (0.48858849999999998, 0.70011546999999996)]),\n",
    "     ('P',\n",
    "      [(6.4648032000000004, 0.15591627),\n",
    "       (1.5022812000000001, 0.60768372000000004),\n",
    "       (0.48858849999999998, 0.39195739000000002)])],\n",
    "    # Ne\n",
    "    [('S',\n",
    "       [(207.01561000000001, 0.15432897000000001),\n",
    "        (37.708151000000001, 0.53532813999999995),\n",
    "        (10.205297, 0.44463454000000002)]),\n",
    "      ('S',\n",
    "       [(8.2463151000000003, -0.099967230000000004),\n",
    "        (1.9162661999999999, 0.39951282999999999),\n",
    "        (0.62322929999999999, 0.70011546999999996)]),\n",
    "      ('P',\n",
    "       [(8.2463151000000003, 0.15591627),\n",
    "        (1.9162661999999999, 0.60768372000000004),\n",
    "            (0.62322929999999999, 0.39195739000000002)])]\n",
    "]\n",
    "basis_set_data = Dict(\"sto3g\" => sto3g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "@assert length(sto3g)==10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Atoms and Molecules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "nat (generic function with 1 method)"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type Atom\n",
    "    atno::Int64\n",
    "    x::Float64\n",
    "    y::Float64\n",
    "    z::Float64\n",
    "end\n",
    "\n",
    "type Molecule\n",
    "    atomlist::Array{Atom,1}\n",
    "end\n",
    "\n",
    "function push!(mol::Molecule,at::Atom)\n",
    "    Base.push!(atomlist,at)\n",
    "end\n",
    "\n",
    "tobohr(x::Float64) = x/0.52918\n",
    "function tobohr!(at::Atom)\n",
    "    at.x /= 0.52918\n",
    "    at.y /= 0.52918\n",
    "    at.z /= 0.52918\n",
    "end\n",
    "function tobohr!(mol::Molecule)\n",
    "    for at in mol.atomlist\n",
    "        tobohr!(at)\n",
    "    end\n",
    "end\n",
    "\n",
    "nuclear_repulsion(a::Atom,b::Atom)= a.atno*b.atno/sqrt(dist2(a.x-b.x,a.y-b.y,a.z-b.z))\n",
    "function nuclear_repulsion(mol::Molecule)\n",
    "    nr = 0\n",
    "    for (i,j) in spairs(nat(mol))\n",
    "        nr += nuclear_repulsion(mol.atomlist[i],mol.atomlist[j])\n",
    "    end\n",
    "    return nr\n",
    "end\n",
    "\n",
    "nel(mol::Molecule) = sum([at.atno for at in mol.atomlist])\n",
    "nat(mol::Molecule) = length(mol.atomlist)\n",
    "\n",
    "# Other molecule methods to implement\n",
    "# nocc, nclosed, nopen, nup, ndown, stoich, mass,\n",
    "# center_of_mass, center!\n",
    "\n",
    "# Array of symbols, masses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample molecules for tests\n",
    "h2 = Molecule([Atom(1,  0.00000000,     0.00000000,     0.36628549),\n",
    "               Atom(1,  0.00000000,     0.00000000,    -0.36628549)])\n",
    "\n",
    "h2o = Molecule([Atom(8,   0.00000000,     0.00000000,     0.04851804),\n",
    "                Atom(1,   0.75300223,     0.00000000,    -0.51923377),\n",
    "                Atom(1,  -0.75300223,     0.00000000,    -0.51923377)])\n",
    "\n",
    "ch4 = Molecule([Atom(6,   0.00000000,     0.00000000,     0.00000000),\n",
    "                Atom(1,   0.62558332,    -0.62558332,     0.62558332),\n",
    "                Atom(1,  -0.62558332,     0.62558332,     0.62558332),\n",
    "                Atom(1,   0.62558332,     0.62558332,    -0.62558332),\n",
    "                Atom(1,  -0.62558332,    -0.62558332,    -0.62558332)])\n",
    "\n",
    "c6h6 = Molecule([ Atom(6,  0.98735329,     0.98735329,     0.00000000),\n",
    "                  Atom(6,  1.34874967,    -0.36139639,     0.00000000),\n",
    "                  Atom(6,  0.36139639,    -1.34874967,     0.00000000),\n",
    "                  Atom(6, -0.98735329,    -0.98735329,     0.00000000),\n",
    "                  Atom(6, -1.34874967,     0.36139639,     0.00000000),\n",
    "                  Atom(6, -0.36139639,     1.34874967,     0.00000000),\n",
    "                  Atom(1,  1.75551741,     1.75551741,     0.00000000),\n",
    "                  Atom(1,  2.39808138,    -0.64256397,     0.00000000),\n",
    "                  Atom(1,  0.64256397,    -2.39808138,     0.00000000),\n",
    "                  Atom(1, -1.75551741,    -1.75551741,     0.00000000),\n",
    "                  Atom(1, -2.39808138,     0.64256397,     0.00000000),\n",
    "                  Atom(1, -0.64256397,     2.39808138,     0.00000000)])\n",
    "\n",
    "lih = Molecule([Atom(3,    0.00000000,     0.00000000,    -0.53999756),\n",
    "                Atom(1,    0.00000000,     0.00000000,     1.08999756)])\n",
    "\n",
    "# Convert to atomic units (bohr)\n",
    "tobohr!(h2)\n",
    "tobohr!(h2o)\n",
    "tobohr!(ch4)\n",
    "tobohr!(c6h6)\n",
    "tobohr!(lih)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "@assert isapprox(nuclear_repulsion(h2),0.7223600367)\n",
    "@assert nel(h2) == 2\n",
    "@assert nel(h2o) == 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dict{Char,Array{Tuple{Int64,Int64,Int64},1}} with 3 entries:\n",
       "  'P' => Tuple{Int64,Int64,Int64}[(1, 0, 0), (0, 1, 0), (0, 0, 1)]\n",
       "  'S' => Tuple{Int64,Int64,Int64}[(0, 0, 0)]\n",
       "  'D' => Tuple{Int64,Int64,Int64}[(2, 0, 0), (0, 2, 0), (0, 0, 2), (1, 1, 0), (…"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type BasisSet # list of CGBFs\n",
    "    bfs::Array{CGBF,1}\n",
    "end\n",
    "\n",
    "basisset() = BasisSet(CGBF[])\n",
    "\n",
    "function push!(basis::BasisSet,cbf::CGBF)\n",
    "    Base.push!(basis.bfs,cbf)\n",
    "end\n",
    "\n",
    "function build_basis(mol::Molecule,name=\"sto3g\")\n",
    "    data = basis_set_data[name]\n",
    "    basis_set = basisset()\n",
    "    for atom in mol.atomlist\n",
    "        for btuple in data[atom.atno]\n",
    "            sym,primlist = btuple\n",
    "            for (I,J,K) in sym2power[sym]\n",
    "                cbf = cgbf(atom.x,atom.y,atom.z,I,J,K)\n",
    "                push!(basis_set,cbf)\n",
    "                for (expn,coef) in primlist\n",
    "                    push!(cbf,expn,coef)\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    return basis_set\n",
    "end\n",
    "\n",
    "sym2power = Dict(\n",
    "    'S' => [(0,0,0)],\n",
    "    'P' => [(1,0,0),(0,1,0),(0,0,1)],\n",
    "    'D' => [(2,0,0),(0,2,0),(0,0,2),(1,1,0),(1,0,1),(0,1,1)]\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "bfs = build_basis(h2)\n",
    "@assert length(bfs.bfs)==2\n",
    "l,r = bfs.bfs\n",
    "@assert isapprox(overlap(l,l),1)\n",
    "@assert isapprox(overlap(r,r),1)\n",
    "@assert isapprox(overlap(l,r),0.6647387449282997)\n",
    "@assert isapprox(kinetic(l,l),0.76003188)\n",
    "@assert isapprox(kinetic(r,r),0.76003188)\n",
    "@assert isapprox(kinetic(l,r),0.24141861181119084)\n",
    "@assert isapprox(coulomb(l,l,l,l), 0.7746059439196398)\n",
    "@assert isapprox(coulomb(r,r,r,r), 0.7746059439196398)\n",
    "@assert isapprox(coulomb(l,l,r,r), 0.5727937653511646)\n",
    "@assert isapprox(coulomb(l,l,l,r), 0.4488373301593464)\n",
    "@assert isapprox(coulomb(l,r,l,r), 0.3025451156654606)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "all_twoe_ints (generic function with 2 methods)"
      ]
     },
     "execution_count": 172,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function all_1e_ints(bfs,mol)\n",
    "    n = length(bfs.bfs)\n",
    "    S = Array{Float64,2}(n,n)\n",
    "    T = Array{Float64,2}(n,n)\n",
    "    V = Array{Float64,2}(n,n)\n",
    "    for (i,j) in pairs(n)\n",
    "        a,b = bfs.bfs[i],bfs.bfs[j]\n",
    "        S[i,j] = S[j,i] = overlap(a,b)\n",
    "        T[i,j] = T[j,i] = kinetic(a,b)\n",
    "        V[i,j] = V[j,i] = nuclear_attraction(a,b,mol)\n",
    "    end\n",
    "    return S,T,V\n",
    "end\n",
    "\n",
    "function all_twoe_ints(bflist,ERI=coulomb)\n",
    "    n = length(bflist.bfs)\n",
    "    totlen = div(n*(n+1)*(n*n+n+2),8)\n",
    "    ints2e = Array{Float64}(totlen)\n",
    "    for (i,j,k,l) in iiterator(n)\n",
    "        #println(\"$i,$j,$k,$l,$(iindex(i,j,k,l))\")\n",
    "        ints2e[iindex(i,j,k,l)] = ERI(bflist.bfs[i],bflist.bfs[j],bflist.bfs[k],bflist.bfs[l])\n",
    "    end\n",
    "    return ints2e\n",
    "end\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dmat (generic function with 1 method)"
      ]
     },
     "execution_count": 173,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function make2JmK(D,Ints)\n",
    "    n = size(D,1)\n",
    "    G = Array{Float64}(n,n)\n",
    "    D1 = reshape(D,n*n)\n",
    "    temp = Array{Float64}(n*n)\n",
    "    for (i,j) in pairs(n)\n",
    "        kl = 1\n",
    "        for (k,l) in rpairs(n)\n",
    "            temp[kl] = 2*Ints[iindex(i,j,k,l)]-Ints[iindex(i,k,j,l)]\n",
    "            kl += 1\n",
    "        end\n",
    "        G[i,j] = G[j,i] = dot(D1,temp)\n",
    "    end\n",
    "    return G\n",
    "end\n",
    "\n",
    "dmat(U,nocc) = U[:,1:nocc]*U[:,1:nocc]'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rhf (generic function with 3 methods)"
      ]
     },
     "execution_count": 174,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function rhf(mol::Molecule,MaxIter::Int64=8,verbose::Bool=false)\n",
    "    bfs = build_basis(mol)\n",
    "    S,T,V = all_1e_ints(bfs,mol)\n",
    "    Ints = all_twoe_ints(bfs)\n",
    "    h = T+V\n",
    "    E,U = eig(h,S)\n",
    "    Enuke = nuclear_repulsion(mol)\n",
    "    nclosed,nopen = divrem(nel(mol),2)\n",
    "    Eold = 0\n",
    "    Energy = 0\n",
    "    println(\"Nel=$(nel(mol)) Nclosed=$nclosed\")\n",
    "    if verbose\n",
    "        println(\"S=\\n$S\")\n",
    "        println(\"h=\\n$h\")\n",
    "        println(\"T=\\n$T\")\n",
    "        println(\"V=\\n$V\")\n",
    "        println(\"E: $E\")\n",
    "        println(\"U: $U\")\n",
    "        println(\"2e ints:\\n$Ints\")\n",
    "    end\n",
    "    for iter in 1:MaxIter\n",
    "        D = dmat(U,nclosed)\n",
    "        if verbose\n",
    "            println(\"D=\\n$D\")\n",
    "        end\n",
    "        G = make2JmK(D,Ints)\n",
    "        H = h+G\n",
    "        E,U = eig(H,S)\n",
    "        Eone = trace2(D,h)\n",
    "        Etwo = trace2(D,H)\n",
    "        Energy = Enuke + Eone + Etwo\n",
    "        println(\"HF: $iter  $Energy : $Enuke    $Eone    $Etwo\")\n",
    "        if isapprox(Energy,Eold)\n",
    "            break\n",
    "        end\n",
    "        Eold  = Energy\n",
    "    end\n",
    "    return Energy,E,U\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Nel=2 Nclosed=1\n",
      "HF: 1  -1.1170995844339022 : 0.7223600367025186    -1.2578170586017612    -0.5816425625346597\n",
      "HF: 2  -1.1170995844339027 : 0.7223600367025186    -1.2578170586017614    -0.5816425625346598\n"
     ]
    }
   ],
   "source": [
    "Energy, E, U = rhf(h2)\n",
    "@assert isapprox(Energy,-1.1170996)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Nel=4 Nclosed=2\n",
      "HF: 1  -7.775292865448999 : 0.9739538361317303    -6.21869911666845    -2.530547584912279\n",
      "HF: 2  -8.016513533651297 : 0.9739538361317303    -6.162738949598296    -2.827728420184731\n",
      "HF: 3  -8.097800662083849 : 0.9739538361317303    -6.108354132748267    -2.9634003654673116\n",
      "HF: 4  -8.136836924932021 : 0.9739538361317303    -6.061698292050672    -3.049092469013079\n",
      "HF: 5  -8.158688062303455 : 0.9739538361317303    -6.023005267164741    -3.109636631270443\n",
      "HF: 6  -8.172474256429597 : 0.9739538361317303    -5.991159038889911    -3.155269053671415\n",
      "HF: 7  -8.181624433516676 : 0.9739538361317303    -5.965138410597584    -3.1904398590508225\n",
      "HF: 8  -8.187678963063593 : 0.9739538361317303    -5.94423397875799    -3.217398820437334\n"
     ]
    },
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mAssertionError: isapprox(Energy, -7.86073270525799)\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mAssertionError: isapprox(Energy, -7.86073270525799)\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:515\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "Energy, E, U = rhf(lih)\n",
    "@assert isapprox(Energy,-7.86073270525799)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HF: 1  34.81226386237448: 9.329511496300162   -2.5990919376852064   28.081844303759524\n",
      "HF: 2  33.618682219475076: 9.329511496300162   -2.414167310513696   26.703338033688606\n",
      "HF: 3  33.569464870941715: 9.329511496300162   -2.3719179753787842   26.61187135002034\n",
      "HF: 4  33.55755807814264: 9.329511496300162   -2.3550117793079934   26.58305836115047\n",
      "HF: 5  33.55393600165691: 9.329511496300162   -2.3469928776750377   26.571417383031786\n",
      "HF: 6  33.55270475466546: 9.329511496300162   -2.3429080743485042   26.566101332713806\n",
      "HF: 7  33.55225631940337: 9.329511496300162   -2.3406996389107078   26.563444462013916\n",
      "HF: 8  33.552085364619266: 9.329511496300162   -2.3394464326558384   26.562020300974943\n",
      "  0.170823 seconds (4.95 M allocations: 124.114 MiB, 7.29% gc time)\n"
     ]
    },
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mAssertionError: isapprox(Energy, -74.9597609118851)\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mAssertionError: isapprox(Energy, -74.9597609118851)\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:515\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "@time Energy,E,U = rhf(h2o)\n",
    "@assert isapprox(Energy,-74.9597609118851)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
